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ABSTRACT 



Context. In the borders of the dead zones of protoplanetary disks, the inflow of gas produces a local density maximum 

that triggers the Rossby wave instability. The vortices that form are efficient in trapping solids. 

Aims. We aim to assess the possibility of gravitational collapse of the solids within the Rossby vortices. 

Methods. We perform global simulations of the dynamics of gas and solids in a low mass non-magnetized self-gravitating 

thin protoplanetary disk with the Pencil Code. We use multiple particle species of radius 1, 10, 30, and 100 cm. The dead 

zone is modelled as a region of low viscosity. 

Results. The Rossby vortices excited in the edges of the dead zone are efficient particle traps. Within 5 orbits after their 
appearance, the solids achieve critical density and undergo gravitational collapse into Mars sized objects. The velocity 
dispersions are of the order of 10 ms^^ for newly formed embryos, later lowering to less than Ims^^ by drag force 
cooling. After 200 orbits, 38 gravitationally bound embryos were formed inside the vortices, half of them being more 
massive than Mars. 



Key words. Keywords should be given 



1. Introduction 

The formation of planets is one of the major unsolved 
problems in modern astrophysics. In the standard core 
accretion scenario, sub-/^(m grains assemble into pro- 
gressively larger bodies through electrostatic interactions 
(Natta et al. 2007), eventually growing into centimeter 
and meter sized boulders. Growth beyond this size, how- 
ever, is halted since these boulders have very poor stick- 
ing properties and are easily destroyed by collisions at the 
velocities assumed to be prevalent in circumstellar disks 
(Benz 2000). Furthermore, centimeter and meter sized 
solids are loosely decoupled from the gas, but remain suf- 
ficiently small to be affected by significant gas drag. The 
resulting headwind from the sub-Keplerian gas reduces 
their angular momentum and forces them into spiral tra- 
jectories onto the star in timescales as short as a few thou- 
sand years (Weidenschilling 1977a). 

A mechanism for overcoming these barriers was pre- 
sented by Kretke & Lin (2007). In the presence of suffi- 
cient ionization, the gaseous disk couples with the am- 
bient weak magnetic field, which triggers the growth of 
the magneto-rotational instability (MRI; Balbus & Hawley 
1991). In its saturated state, a vigorous turbulence drives 
accretion onto the star by means of magnetic and kinetic 
stresses. However, in the water condensation front (snow- 
line) the abimdant presence of snowflakes effectively re- 
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moves free electrons from the gas, lowering the degree of 
ionization. The turbulence is weakened locally and the ac- 
cretion flow is stalled. As the radial inflow proceeds from 
the outer disk, gas accumulates at the snowline. Since em- 
bedded solid bodies move towards gas pressure maxima 
(Haghighipour & Boss 2003), the snowline environment 
proposed by Kretke & Lin (2007) is potentially an efficient 
particle trap. This scenario was further explored by Brauer 
et al. (2008), who demonstrated that as solids concentrate 
at this local pressure maximum, rapid growth into kilome- 
ter sized planetesimals occurs by coagulation. 

Kretke et al. (2008) emphasized that an identical 
mechanism is supposed to occur elsewhere in the disk. 
Ionization ought to be present in the very inner disk due 
to the high temperatures, as well as in the outer regions 
where the gas is sufficiently thin for cosmic rays to pene- 
trate to the disk midplane and provide ionization through- 
out. In between, however, temperatures are too low and 
gas columns too thick to allow sufficient ionization either 
by collisions or by cosmic rays. In the midplane of this re- 
gion, the gas is neutral and the turbulence is largely sup- 
pressed (Gammie 1996). As in the snowline, the accretion 
flow from the MRI-active regions halts at the borders of 
this "dead" zone, where the gas then accumulates. 

These models have been tested only by one- 
dimensional simulations, and these tests have therefore 
not benefitted from an interesting development. As 
shown by Varniere & Tagger (2006), the density pileup 
at the border of the dead zone triggers the Rossby wave 
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instability (RWI; Li et al. 2001). The azimuthal symmetry 
of the problem is broken and long-lived anticyclonic 
vortices are formed as the waves break and coalesce. Such 
entities are of significant interest because, by rotating 
clockwise in the global counterclockwise Keplerian flow, 
they amplify the local shear and induce a net force on 
solid particles towards their center (Barge & Sommeria 
1995). As shown by Klahr & Bodenheimer (2006), the 
accumulation of solids under these circumstances is likely 
to lead to high densities. Inaba & Barge (2006) contin- 
ued the study of Varniere & Tagger (2006) by including 
solids and confirming that the Rossby vortices excited 
in the borders of the dead zone act as powerful traps, 
enhancing the local solids-to-gas ratio by at least an order 
of magnitude. Unfortunately, they used a fluid approach 
- which limited the maximum particle size they could 
consider - and they did not include the self-gravity of the 
solids, which is crucial to follow the gravitational collapse. 
Studies with interacting particles in the literature include 
a MRI-unstable local box (Johansen et al. 2007) capable 
of producing dwarf planets out of meter sized boulders, 
a global massive disk imstable to gas self-gravity (Rice 
et al. 2006) in which concentrations of 0.5 M® are seen 
in the spiral arms, and a simulation that produces 10-100 
km "sandpile" clumps formed out of mm-sized particles 
(Cuzzi et al. 2008). 

In this Letter, we build on the studies of Varniere & 
Tagger (2006) and Inaba & Barge (2006) by including self- 
gravitating centimeter and meter sized Lagrangian parti- 
cles to model the solid phase. We show that in the vortices 
launched by the RWI in the borders of the dead zone, the 
solids quickly achieve critical densities and undergo grav- 
itational collapse into protoplanetary embryos in the mass 
range 0.1-0.6 M®. 

2. Model 

We work in the thin disk approximation, using the verti- 
cally averaged equations of hydrodynamics. The gas drag 
is implemented in the same way as Paardekooper (2007), 
interpolating between Epstein and Stokes drag (see on- 
line supplement). The back reaction of the drag force onto 
the gas is present. The Poisson solver for the particles is a 
particle-mesh solver based on multiple Fourier transforms 
in a Cartesian grid, as used by Johansen et al. (2007). 

We follow the Varniere & Tagger (2006) dead zone 
model, which consists of jumps in the viscosity profile. 
We artificially place the inner and outer edges of the dead 
zone at 0.6 and 1.2 times the semi-major axis of Jupiter 
(5.2 AU), using Heaviside functions to jump from a = 
10~^ to zero inside the dead zone. The parameter a is the 
usual alpha viscosity (Shakura & Sunyaev 1973). 

The simulations were performed with the PENCIL 
CODE^ in Cartesian geometry. The resolution was 
256 X 256. The surface density profile followed a power law 
of I^{r) (X r^^/^, and we chose a disk about twice as mas- 
sive as the Minimum Mass Solar Nebula, with surface den- 
sity 1^0=300 gcm^^ at 5.2 AU. The sound speed followed 
the local isothermal approximation with a radial tempera- 
ture profile r(r) ex r^^. The disk aspect ratio was /z=0.05. 
For the solids, we used 10^ Lagrangian superparticles and 

^ See http: / / www.nordita.org/ software / pencil-code 



the interstellar value for the solids-to-gas ratio (10 ^). Each 
superparticle therefore contained IO^^Mq ~ 2 x 10^^ g of 
material. We used multiple particles species, of 1, 10, 30, 
and 100 cm radii, each represented by 1 /4 of the total num- 
ber of particles. We quote time as orbital periods at 5.2 AU. 

3. Results 

In Fig. 1, we show snapshots of the appearance of the disk 
for the gas and solids phases. The vortices triggered by 
the Rossby wave instability are visible as early as 15 or- 
bits. As seen in the solids phase, the particles are trapped 
by the vortical motion and soon reach extremely high den- 
sities. After 17 orbits, seven vortices appear at the outer 
edge. After 45 orbits have elapsed, the m=4 mode begins 
to dominate, persisting until the end of the simulation at 
200 orbits, their gas surface density peaking at 4.5 times 
their initial value. In the inner edge of the dead zone, at 40 
orbits we see a conspicuous m=3 mode. By the end of the 
simulation, their surface density has increased by a factor 
of 8 relative to the initial condition, and a weak m=2 mode 
is visible, albeit with far less contrast than in the outer disk. 

In Fig. 2a, we plot the time evolution of the maxi- 
mum concentration of solids. The solid line represents 
the maximum mass of solids contained in a single cell. 
The red dashed line marks the maximum mass that is 
gravitationally bound. We decide whether boundness is 
present based on two criteria. First, we consider the mass 
inside the Hill's sphere of the clump defined by the 
black line. Particles inside /outside the Hill's sphere are 
added/removed from the total mass, and the Hill's radius 
recomputed. The process is iterated until convergence. 
This positional criterion is followed by a dynamical one. 
We calculate the velocity dispersion Vrms of the particles 
inside the Hill's radius, and compare its value with the es- 
cape velocity of the enclosed mass. If i'rms< ^esc/ we con- 
sider that the cluster of particles is bound. We plot the ve- 
locity dispersion and escape velocity of the most massive 
clump in Fig. 2b. The first bound clumps appeared at 18 or- 
bits, with masses of 0.050 and 0.036 M®. At 20 orbits, four 
clumps of 0.1 M® are present. The mass is that of Mars, 
constituting a protoplanetary embryo. The efficiency of the 
vortex trapping mechanism can be more clearly appreci- 
ated if we consider the time elapsed between the rise of 
the Rossby vortices and the collapse of the trapped parti- 
cles into a Mars sized object: only 5 orbits. 

Two orbits later, the maximum mass increased to 0.22 
M(jj. Nine other clumps collapsed into embryos as well, 
five of them of mass above that of Mars. The maximum 
bound mass reached 0.67 M®, but settled at 0.59 M® from 
89 orbits until the end of the simulation at 200 orbits. We 
observe evidence that the mass loss episodes are due to 
tides from the gas, since the vortices concentrate sufficient 
gas to provide a considerable gravitational pull. In addi- 
tion to tides, erosion (Cuzzi et al. 2008) also plays a role 
in disrupting clumps of smaller particles. We discuss this 
point further in the online supplement. 

Figure 2b also indicates that the velocity dispersion re- 
mained below 1 ms^^ for most of the simulation. This 
is of extreme importance because it implies that particle 
encounters are gentle enough for destructive collisions to 
be avoided. The opposite was reported for the massive 
disk models of Rice et al. (2006), where particle encoun- 
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Fig. 1. The appearance of the disk in the gas (upper panels) and solid (lower panels) phases in selected snapshots. The Rossby vortices 
first appear at 15 orbits. Collapse of the particles into a gravitationally bound planetary embryo the size of Mars occurs 5 orbits later. 
The vortices tend to merge and decrease in number, retaining an m=4 dominant mode until the end of the simulation, up to which 
over a dozen embryos were formed. 
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Fig. 2. a). Time evolution of the maximum mass of solids. The black solid line is the maximum mass of solids inside a single cell, and 
the red dot-dashed line marks the mass that is gravitationally bound. The maximum mass settles at 0.59 Mg^ from 90 orbits onwards. 

b) . The internal velocity dispersion Dj-ms of the most massive embryo, compared with its escaped velocity fesc- The embryo is so 
firmly bound that even the maximum internal speed Umax is lower than the escape velocity. Throughout most of the simulation, the 
velocity dispersion is below 1 m s^^ . 

c) . The mass spectrum by the end of the simulation. The line and color style is the same as in Fig. 2a. Over a dozen gravitationally 
bound embryos in the mass range 0.1-0.6 M® were formed in the vortices launched at the edges of the dead zone. 



ters in the spiral arms occurred at velocities comparable 
with the sound speed. In Fig. 2b, we plot the maximum 
speed fmax for comparison. It is evident that even Umax re- 
mains lower than the escape velocity. This indicates that 
even if destructive collisions occur, the fragments will re- 
main bound, although the strong drag force felt by the 
fragments might delay any gravitational collapse. 

We plot the mass spectrum of the formed embryos in 
Fig. 2c. The solid black line represents the mass of solids 
inside a cell of the simulation box, without information on 
boimdness. We overplot the distribution of bound clumps 
with the red dashed line. Thirty-eight bound embryos are 
formed by the end of the simulation, twenty of these being 
more massive than Mars. The two most massive embryos 
have 0.59 and 0.57 M^, and are located in the inner and 
outer edge, respectively. 



A common trait of these embryos is that they consist of 
an overwhelming majority of particles of a single species. 
For instance, over 98% of the most massive embryo con- 
sists of 30 cm particles. Almost 97% of the second most 
massive one were 10 cm particles. In one embryo in the 
Mars mass bin, 94% of the particles were of 1 cm radius. 
These different particle radii reflect the variations in gas 
drag strengths experienced by particles of different size. 
Since particles of a given size drift in a similar way, their 
spatial convergence becomes easier than the convergence 
of particles of different sizes. The almost single-phasing 
of the size spectrum of the bound clumps also contributes 
to their low velocity dispersions. Since particles of differ- 
ent sizes experience different gas drags, their relative ve- 
locities are often far higher than the relative velocities be- 
tween particles of identical size (Thebault et al. 2008). We 
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Fig. 3. Gravitational collapse of the particles into an embryo twice 
the mass of Mars. The unit of length is the Earth-Moon mean 
separation (-Rem)- The lower left panel has the grid scale over- 
plotted. The circle marks the Hill's radius. A control run without 
selfgravity attains higher cell concentration due to the absence 
of gas tides. The bar shows the mass of the embryo in units of 
M]v[ars ^rid the fractional mass of each particle species. 

indeed see that clumps with a more non-uniform distribu- 
tion of particle sizes usually show larger velocity disper- 
sions (~10 instead of ~1 ms^^). It is an observed fact that 
planetesimals are formed of similarly sized particles (Scott 
& Krot 2005). Although the constituent particles appear to 
be sub-mm grains, different nebula parameters could en- 
able smaller grains to be trapped. 

At higher resolution (512^, and 4 x 10^ particles), more 
bound clumps are formed, which extend to lower masses. 
More clumps containing similar quantities of particles 
of different sizes are observed, although the majority of 
clumps have nearly single phasing. The situation does 
not alter significantly when the number of particles is in- 
creased to 10^. In Fig. 3, we follow the collapse of one 
of these clumps at higher resolution. Although particles 
of different size are present inside the Hill's sphere dur- 
ing the first stages of collapse, most are expelled and the 
collapsed embryo consists primarily of 30 cm particles. A 
control run without self-gravity achieves higher cell con- 
centration, due to the absence of gas tides. At later times, 
the clumps also tend to have low rms speeds and small 
rms radii (falO"* km), due to the efficient dynamical cool- 
ing provided by the drag force. Size segregation due to 
aerodynamical sorting is also observed in the control run. 
Inaba & Barge (2006) reported destruction of the vortices 
by the drag force backreaction. We see a different effect, in 
which the particles alter the evolution of the Rossby vor- 
tices and generate a vorticity of their own. 

4. Conclusions 

In this Letter, we have shown that when modelling the 
self-gravity of the gas and solids in protoplanetary disks, 
gravitational collapse of the solids into Mars sized pro- 
toplanetary embryos occurs rapidly at the borders of the 
dead zone, where particles concentrate. We have also 
found that tides from the dense gaseous vortices may hin- 
der the formation process significantly 



Studies considering the origin of oligarchs usually be- 
gin from a collection of 10-20 Mars sized objects (e.g., 
Kobuko et al. 2006). This Letter presents the first simu- 
lation in which a similar number of Mars-sized embryos 
are formed from centimeter and meter sized sized build- 
ing blocks. 

It is crucial to the model that particles grow to suffi- 
cient size, otherwise the drag force from the gas becomes 
too strong to allow any concentration. Testi et al. (2003) 
observed grains of up to cm sizes in the disk surround- 
ing the pre-main-sequence star CQ Tauri, which provides 
some observational evidence that a sufficient number of 
particles of the required size may exist in true protoplane- 
tary disks. 

We emphasize again that the model used in this Letter 
is simplistic, and the conditions may differ with a more re- 
alistic treatment of the dead zone. Nevertheless, the mech- 
anism presented in this Letter (as proposed originally by 
Varniere & Tagger 2006) appears robust. It only requires 
the RWI to be excited in the borders between the active and 
dead zones, which in turn relies only on the slowdown 
of the accretion flow at this same border. Future studies 
should include a coagulation/fragmentation model such 
as those of Brauer et al. (2008) or Johansen et al. (2008), and 
focus on the precise state of flow at this transition region 
in global simulations, to address the question of how the 
RWI interacts with the MRI and non-ideal MHD in three 
dimensions. 
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Appendix A: Dynamical Equations 

As stated in the main paper, we work in the thin disk ap- 
proximation, using the vertically integrated equations of 
hydrodynamics 



dt 

du 
dt 
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In the above equations, G is the gravitational constant. 
Eg and Ep are the vertically integrated gas density and 
bulk density of solids, respectively, u represents the ve- 
locity of the gas parcels, Xp is the position and Vp is the 
velocity of the solid particles, P is the vertically integrated 
pressure, Cg is the sound speed, <P is the gravitational po- 
tential, V is the viscosity, and S is the rate of strain ten- 
sor. The functions fi){E^) and fv{u, E^) are sixth order hy- 
perdiffusion and hyperviscosity terms that provide extra 
dissipation close to the grid scale, explained in Lyra et al. 
(2008). They are required because the high order scheme of 
the Pencil code has too little overall numerical dissipation. 

The ftmction is the drag force by which gas and 
solids interact. In Eq. (A.8), p, is the internal density of a 
solid particle, fl, its radius, and Av = Vp — u its velocity 
relative to the gas. Cp is a dimensionless coefficient that 
defines the strength of the drag force. 



Appendix B: Drag Force 

The drag regimes are controlled by the mean free path A of 
the gas, which can be expressed in terms of the Knudsen 
number of the flow past the particle Kn = A/2fl,. High 
Knudsen numbers correspond to free molecular flow, or 
Epstein regime. Stokes drag is applicable to low Knudsen 
numbers. We use the formula of Woitke & Helling (2003; 
see also Paardekooper 2007), which interpolates between 
Epstein and Stokes regimes 



Cr 



where and Cp* are the coefficients of Epstein and 
Stokes drag, respectively. These coefficients are 
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where Ma = |Ap|/cs is the Mach number. Re — 
2a,p\Av\/}i is the Reynolds number of the flow past the 
particle, and }i = y/STrrpCsA/S is the kinematic viscosity 
of the gas. 

The approximation for Epstein drag (Kwok 1975) 
connects regimes of low and high Mach number with a 
good degree of accuracy, and is more numerically friendly 
than the general case (Baines et al. 1965). The piecewise 
function for the Stokes regime are empirical corrections 
to Stokes law (Cd = 24Re-^), which only applies for low 
Reynolds numbers. 



Appendix C: Gas tides and mass loss 

The most remarkable feature of Fig. 3a of the main paper 
is the mass loss event at 90 orbits. It consists of the de- 
tachment of a 0.8 MMars object from the original cluster, of 

6.7MMars- 

We observe evidence that this puzzling behaviour is 
due to gravitational tides from the gas. The gas is too 
pressure-supported to undergo collapse, but the vortices 
concentrate a sufficient amount of material to exert a con- 
siderable gravitational pull. We illustrate this in Fig. C.l, 
where we show the state of the disk before (at 80 orbits. 
Fig. C. la-Fig. C.ld) the mass-loss episode and after (at 
100 orbits. Fig. C.le-Fig. C.lh). The plots show the sur- 
face densities of gas and solids, and the potential associ- 
ated with them. Even though the clumping of solids yield 
a considerable gravitational pull (Fig. C.ld and Fig. C.lh), 
these figures indicate that the dominant contribution to the 
gravitational potential comes from the gas - more specifi- 
cally from the vortices, where the gas density peaks at a 
value one order of magnitude higher than that of the ini- 
tial conditions. 

The most massive clump is located in the inner disk 
at (X,Y)=(-0.40,-0.53) in Fig. C.lb, not clearly identifi- 
able amidst the other particles trapped inside the vortex. 
However, the embryo is immediately observable as the 
bright point at (X,Y)=(-0.65,-0.19) in Fig. C.lh (also visi- 
ble in Fig. C.lf, albeit less prominently). There are two fea- 
tures in this plot that are worth noting. First, by compar- 
ing the location of the embryo in these figures with the 
location of the vortices, we notice that the planet has left 
its parental vortex. Second, the inner vortices have under- 
gone the transition from the m=3 to the m=2 mode. Due 
to merging, their gas density has increased, with dramatic 
consequences for the embryos within them. 

We assess how the gravity of the gas influences the mo- 
tion of the particles in Fig. C.li-Fig. C.ll. In Fig. C.li, we 
take a horizontal density slice at the position of the most 
massive embryo at 80 orbits. Figure C.lj is an enlargement 
of Fig. C.li about X=-0.53, where the embryo is located. 
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Fig. C.l. The state of the disk before {a-d) and after {e-h) the mass-loss episode. The conspicuous differences between them are due 
to the inner vortices passing from the m=3 to the m=2 mode. As a result, the the most massive embryo has left its parental vortex 
between one snapshot and another. It is seen as a bright spot in panels / and h, at (X,Y)=(-0. 65,-0. 19). In panel h (before the mass loss), 
the embryo is at (X,Y)=(-0.40,-0.53) but not easily identified among the swarm of solids inside the vortex. Panel i shows a horizontal 
slice through this location, in which we see that the peak in the density of solids is not much higher than the gas density at the location 
of the embryo (panel j). Significant gas tides are expected because the gravitational potential (panel k) and acceleration (panel /) have 
similar contributions from gas and solid components. 



We see that the densities of solids and gas peak at simi- 
lar values. The subsequent figures show the gravitational 
potential (Fig. C.Ik) and acceleration (Fig. C.ll) about the 
embryo. The gas produces a deeper gravitational well, al- 
beit smoother than the one displayed by the solids. In the 
acceleration plot, it is seen that the pull of the gas is more 



significant than the pull of the embryo at a distance of only 
0.26 AU (0.03 in code units, corresponding to two grid 
cells) away from the center. And even where the pull of 
the solids is strongest (one grid cell away from the center), 
the gravity of the gas is still an appreciable fraction of the 
gravity of the solids. Tides from the gas are imavoidable. 
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It is beyond the scope of this paper to consider the full 
mathematical details of the theory of tides, especially be- 
cause the two bodies (the vortex and the embryo) are ex- 
tended. Instead, we consider the following toy model. The 
tidal force Fj experienced by the planet is proportional to 
the gradient of the acceleration a induced by the gas. It is 
also proportional to the radius R of the planet: Fj c< R Vfl. 
Since Va = —V^<P, according to the Poisson equation, the 
tidal force is proportional to the local value of the density 



FjocRpg. (C.l) 

We consider the 3D volume density to avoid the 
requirement of using the Dirac delta in the 2D case. 
Considering the planet spherical, Newton's second the- 
orem holds and we can assume that Fc = -GM/R^ 
for the planet's (self-)gravitational force at its surface. 
Substituting M = i/3nppR^, we have Fq cx Rpp, so 

^ = ^ - ^, (C.2) 

Pp 

i.e., the ratio of the disrupting tidal stresses to the self- 
gravitating forces that attempt to keep the planet together 
is directly proportional to the gas-to-solids ratio. At 80 or- 
bits, as seen in Fig. C.lj, this ratio is around unity. As the 
vortices undergo the transition from the m=3 to the m=2 
mode, their peak density increases (while the planet re- 
mains at constant mass). The tides eventually become suf- 
ficiently strong for a mass-loss event to occur. We also wit- 
ness some of the other, less massive, embryos being com- 
pletely obliterated by the gas tides. Erosion is also impor- 
tant, since we observe a size dependency in this effect, 
with embryos consisting of a, =10 cm particles being more 
prone to disruption than those consisting of fl,=30 cm par- 
ticles. We performed tests that indicated that the erosion of 
bound clumps by ram pressure seen by Cuzzi et al. (2008) 
only occurs for clumps consisting of particles smaller than 
cm-size for our initial nebula parameters, larger particles 
being more stable. However, when the gas density of the 
vortices reaches a maximum value a factor ten higher than 
the initial density, the gas drag also strengthens, shifting 
the onset of erosion towards larger particle radii. 

The effect of tides will probably be less dramatic in 3D 
simulations because, as the particles settle into the mid- 
plane, the ratio of the volume gas density to the bulk den- 
sity of solids p<^/pp is expected to be much lower than the 
ratio of the column gas density to the vertically integrated 
surface density of solids Eg/Lp. 



